LANDSAT

Freely available: https://earthexplorer.usgs.gov/

Different bands

Load imagery and create raster brick

library(raster)
library(sp)

l.blue <- raster('LC08_L1TP_024030_20170602_20170615_01_T1/LC08_L1TP_024030_20170602_20170615_01_T1_B2.TIF')
l.green <- raster('LC08_L1TP_024030_20170602_20170615_01_T1/LC08_L1TP_024030_20170602_20170615_01_T1_B3.TIF')
l.red <- raster('LC08_L1TP_024030_20170602_20170615_01_T1/LC08_L1TP_024030_20170602_20170615_01_T1_B4.TIF')
l.nir <- raster('LC08_L1TP_024030_20170602_20170615_01_T1/LC08_L1TP_024030_20170602_20170615_01_T1_B5.TIF')
l.swir1 <- raster('LC08_L1TP_024030_20170602_20170615_01_T1/LC08_L1TP_024030_20170602_20170615_01_T1_B6.TIF')
l.swir2 <- raster('LC08_L1TP_024030_20170602_20170615_01_T1/LC08_L1TP_024030_20170602_20170615_01_T1_B7.TIF')
l.pc <- raster('LC08_L1TP_024030_20170602_20170615_01_T1/LC08_L1TP_024030_20170602_20170615_01_T1_B8.TIF')
l.cirrus <- raster('LC08_L1TP_024030_20170602_20170615_01_T1/LC08_L1TP_024030_20170602_20170615_01_T1_B9.TIF')
l.t1 <- raster('LC08_L1TP_024030_20170602_20170615_01_T1/LC08_L1TP_024030_20170602_20170615_01_T1_B10.TIF')
l.t2 <- raster('LC08_L1TP_024030_20170602_20170615_01_T1/LC08_L1TP_024030_20170602_20170615_01_T1_B11.TIF')

# Create raster brick
ls = brick(l.red,l.green,l.blue,l.nir,l.swir1,l.swir2,l.cirrus,l.t1,l.t2)

# CROP
e <- extent(280000, 330000, 4750000, 4800000)
ls.c = crop(ls,e)

names(ls.c) <- c('red','green','blue','NIR','SWIR1','SWIR2','cirrus','thermal1','thermal2')
writeRaster(ls.c,filename = 'Madison_Landsat_20170602.tif', format="GTiff", overwrite=TRUE)

Load cropped imagery

library(raster)
## Loading required package: sp
library(sp)
r1 = brick('Madison_Landsat_20170602.tif')
r2 = brick('Madison_Landsat_20180317.tif')
names(r1) <- c('red','green','blue','NIR','SWIR1','SWIR2','cirrus','thermal1','thermal2')
names(r2) <- c('red','green','blue','NIR','SWIR1','SWIR2','cirrus','thermal1','thermal2')

r1
## class       : RasterBrick 
## dimensions  : 1667, 1667, 2778889, 9  (nrow, ncol, ncell, nlayers)
## resolution  : 30, 30  (x, y)
## extent      : 280005, 330015, 4750005, 4800015  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=16 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : /Users/hilarydugan/Documents/Zoo955/Lecture9_Landsat/Madison_Landsat_20170602.tif 
## names       :   red, green,  blue,   NIR, SWIR1, SWIR2, cirrus, thermal1, thermal2 
## min values  :  3868,  5598,  7421,  5232,  5142,  5058,   4995,    20368,    19692 
## max values  : 65535, 65535, 61486, 65535, 65535, 65535,   6626,    38316,    32482

Plot imagery

plot(r2)

Plot RGB

plotRGB(r1, r = 1, g = 2, b = 3, axes = TRUE, stretch = "lin",
        main = "Landsat True Color Composite")

plotRGB(r2, r = 1, g = 2, b = 3, axes = TRUE, stretch = "lin",
        main = "Landsat True Color Composite")

Extract values

# Load Lake shapefile
library(sf)
## Warning: package 'sf' was built under R version 3.4.4
## Linking to GEOS 3.6.1, GDAL 2.1.3, proj.4 4.9.3
mendota <- st_read('../Lecture2_CRS/Data/LakeMendota.shp') %>%
  st_transform('+proj=utm +zone=16 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0')
## Reading layer `LakeMendota' from data source `/Users/hilarydugan/Documents/Zoo955/Lecture2_CRS/Data/LakeMendota.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 12 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -89.48369 ymin: 43.07384 xmax: -89.36737 ymax: 43.14706
## epsg (SRID):    4269
## proj4string:    +proj=longlat +datum=NAD83 +no_defs
# random sample of 100 points
l.pts <- st_sample(mendota, size = 100) %>% st_sf
lake.s <- extract(r1, as(l.pts,'Spatial')) #summer values
lake.w <- extract(r2, as(l.pts,'Spatial')) #winter values

# Load watershed shapefile
watershed <- st_read('../Lecture3_Shapefiles/Data/YaharaBasins/Mendota_Basin.shp') %>%
  st_transform('+proj=utm +zone=16 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0')
## Reading layer `Mendota_Basin' from data source `/Users/hilarydugan/Documents/Zoo955/Lecture3_Shapefiles/Data/YaharaBasins/Mendota_Basin.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 3 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 552122.5 ymin: 286231.1 xmax: 584702.5 ymax: 321699.8
## epsg (SRID):    NA
## proj4string:    +proj=tmerc +lat_0=0 +lon_0=-90 +k=0.9996 +x_0=520000 +y_0=-4480000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
# random sample of 100 points
ws.pts <- st_sample(watershed, size = 100) %>% st_sf
ws.s <- extract(r1, as(ws.pts,'Spatial'))
ws.w <- extract(r2, as(ws.pts,'Spatial'))

# To see some of the reflectance values
lake.S <- colMeans(lake.s)
lake.W <- colMeans(lake.w)
ws.S <- colMeans(ws.s,na.rm = T)
ws.W <- colMeans(ws.w,na.rm = T)

plot(ws.S,xlab='Bands',ylab='Reflectance',type='b',xaxt='n')
axis(1,at = 1:9,labels = names(r2),las=2)

lines(lake.S,type='b',col='blue')
lines(ws.W,type='b',col='black',pch=16)
lines(lake.W,type='b',col='blue',pch=16)
legend('topleft',legend = c('Watershed - Summer','Lake - Summer','Watershed - Winter','Lake - Winter'),col = c('black','blue','black','blue'),
       pch=c(21,21,16,16),bty='n')